home *** CD-ROM | disk | FTP | other *** search
/ Mac Power 1997 December / MACPOWER-1997-12.ISO.7z / MACPOWER-1997-12.ISO / AMUG / PROGRAMMING / Raven 1.2.sit / Raven 1.2 / Source / Foundation / Common / ZNumbers.cpp < prev    next >
Text File  |  1997-06-21  |  13KB  |  511 lines
  1. /*
  2.  *  File:       ZNumbers.cpp
  3.  *  Summary:       Number related functions.
  4.  *  Written by: Jesse Jones
  5.  *
  6.  *  Copyright ｩ 1996-1997 Jesse Jones. 
  7.  *    For conditions of distribution and use, see copyright notice in ZTypes.h  
  8.  *
  9.  *  Change History (most recent first):    
  10.  *
  11.  *         <4>     4/27/97    JDJ        Added int and double versions of Random.
  12.  *         <3>     4/15/97    JDJ        Renamed EqualReal Equal.
  13.  *         <2>     2/23/97    JDJ        EqualReal ASSERTs that tolerance and typical are valid.
  14.  *         <1>     1/13/96    JDJ        Created
  15.  */
  16.  
  17. #include <ZNumbers.h>
  18.  
  19. #include <Events.h>
  20. #include <Fp.h>
  21. #include <Limits.h>
  22.  
  23. #include <ZDebug.h>
  24. #include <ZTypes.h>
  25. #include <ZUnitTest.h>
  26.  
  27.  
  28. // ===================================================================================
  29. //    Internal Functions
  30. //         The random number algorithm  is from D. E. Knuth,  _The Stanford GraphBase_, 
  31. //        Addison-Wesley, ISBN 0-201-54275-7. Note that the generator is much better 
  32. //        than a simple linear congruential algorithm, but about 13% slower.
  33. // ===================================================================================
  34.  
  35. long    gb_flip_cycle(void);
  36. long    gb_init_rand(long seed);
  37. long    gb_unif_rand(long m);
  38. int        gb_flip(void);
  39.  
  40. const unsigned long    kTwo_to_the_31 = 0x80000000UL;
  41.  
  42.  
  43. //---------------------------------------------------------------
  44. //
  45. // Private Declarations
  46. //
  47. //---------------------------------------------------------------
  48. static long A[56] = {-1};
  49.  
  50. long sFlipSeed = gb_init_rand((long) TickCount());        
  51.  
  52.  
  53. //---------------------------------------------------------------
  54. //
  55. // External Declarations
  56. //
  57. //---------------------------------------------------------------
  58. long *gb_fptr = A;
  59.  
  60. #define gb_next_rand() (*gb_fptr >= 0 ? *gb_fptr-- : gb_flip_cycle())
  61.  
  62.  
  63. #define mod_diff(x,y) (((x) - (y)) & 0x7FFFFFFF)
  64.  
  65.  
  66. //---------------------------------------------------------------
  67. //
  68. // gb_flip_cycle
  69. //
  70. //---------------------------------------------------------------
  71. long gb_flip_cycle(void)
  72. {
  73.     register long *ii, *jj;
  74.     
  75.     for (ii = &A[1], jj = &A[32]; jj <= &A[55]; ii++, jj++)
  76.            *ii = mod_diff(*ii, *jj);
  77.     for (jj = &A[1]; ii<= &A[55]; ii++, jj++)
  78.            *ii = mod_diff(*ii, *jj);
  79.     gb_fptr = &A[54];
  80.     return A[55];
  81. }
  82.  
  83.  
  84. //---------------------------------------------------------------
  85. //
  86. // gb_init_rand
  87. //
  88. //---------------------------------------------------------------
  89. long gb_init_rand(long seed)
  90. {
  91.     register long i;
  92.     register long prev = seed, next = 1;
  93.     
  94.     seed = prev = mod_diff(prev, 0);
  95.     A[55] = prev;
  96.     for (i = 21; i; i = (i + 21) % 55) {
  97.         A[i] = next;
  98.         /* Compute a new next value. */
  99.         next = mod_diff(prev, next);
  100.         if (seed & 1)
  101.             seed = 0x40000000 + (seed >> 1);
  102.         else
  103.             seed >>= 1;
  104.         next = mod_diff(next, seed);
  105.  
  106.         prev = A[i];
  107.     }
  108.     /* Get the array values ``warmed up''. */
  109.     (void) gb_flip_cycle();
  110.     (void) gb_flip_cycle();
  111.     (void) gb_flip_cycle();
  112.     (void) gb_flip_cycle();
  113.     (void) gb_flip_cycle();
  114.     
  115.     if (seed == 0)
  116.         seed = 1;
  117.     
  118.     return seed;                                    // Hack to enable gb_init_rand to be called before main.
  119. }
  120.  
  121.  
  122. //---------------------------------------------------------------
  123. //
  124. // gb_unif_rand
  125. //
  126. //---------------------------------------------------------------
  127. long gb_unif_rand(long m)
  128. {
  129.     ASSERT(m > 0);            
  130.     
  131.     register long r;
  132.     register unsigned long t = kTwo_to_the_31 - (kTwo_to_the_31 % m);
  133.     
  134.     if (m > 1) {
  135.         do {
  136.             r = gb_next_rand();
  137.         } while (t <= (unsigned long) r);
  138.         r = r % m;
  139.         
  140.     } else
  141.         r = 0;
  142.     
  143.     return r;
  144. }
  145.  
  146.  
  147. //---------------------------------------------------------------
  148. //
  149. // gb_flip
  150. //
  151. //---------------------------------------------------------------
  152. int gb_flip(void) 
  153. {
  154.     int j;
  155.     
  156.     gb_init_rand(-314159L);
  157.     if (gb_next_rand() != 119318998) {
  158.         DEBUGSTR("Failure on the first try.¥n");
  159.         return -1;
  160.     }
  161.     for (j = 1; j <= 133; j++) 
  162.         gb_next_rand();
  163.     if (gb_unif_rand(0x55555555L) != 748103812) {
  164.         DEBUGSTR("Failure on the second try.¥n");
  165.         return -2;
  166.     }
  167.     /* DEBUGSTR("OK, the gb_flip_routines seem to work.¥n"); */
  168.     return 0;
  169. }
  170.  
  171. #pragma mark -
  172.  
  173. // ===================================================================================
  174. //    Global Functions
  175. // ===================================================================================
  176.  
  177. //---------------------------------------------------------------
  178. //
  179. // FlipCoin
  180. //
  181. // The algorithm is from the second edition of _Numerical Recipies 
  182. // in C_ by Press, Teukolsky, Vetterling, and Flannery and is based
  183. // on "primitive polynomials modulo 2". The polynomial used here
  184. // is x**18 + x**5 + x**2 + x + 1. This is guaranteed to produce 
  185. // all possible sequences of 18 bits before repeating.
  186. //
  187. // Numerical Recipies has a table with larger polynomials that
  188. // will yield a larger period but they'll also yield (occasionally)
  189. // longer runs returning the same result.
  190. //
  191. //---------------------------------------------------------------
  192. bool FlipCoin()
  193. {
  194.     const ulong kBit1  = 1;
  195.     const ulong kBit2  = 2;
  196.     const ulong kBit5  = 16;
  197.     const ulong kBit18 = 131072;
  198.     
  199.     const ulong kMask = kBit1 + kBit2 + kBit5;
  200.     
  201.     ASSERT(sFlipSeed != 0);
  202.     
  203.     if (sFlipSeed & kBit18) {
  204.         sFlipSeed = (long) (((sFlipSeed ^ kMask) << 1) | kBit1);
  205.         return true;
  206.         
  207.     } else {
  208.         sFlipSeed <<= 1;
  209.         if (sFlipSeed == 0)            
  210.             sFlipSeed = 1;
  211.         return false;
  212.     }
  213. }
  214.  
  215.  
  216. //---------------------------------------------------------------
  217. //
  218. // Equal
  219. //
  220. //---------------------------------------------------------------
  221. bool Equal(double x, double y, double tolerance, double typical)
  222. {
  223.     ASSERT(tolerance >= 0.0);
  224.     ASSERT(typical > 0.0);
  225.     
  226.     bool equal = false;
  227.     
  228.     double z = Max(Abs(x), Abs(y));
  229.     
  230.     if (z == __inf())
  231.         equal = x == y;
  232.  
  233.     if (!equal) {
  234.         if (Abs(x - y)/z <= tolerance)
  235.             equal = true;
  236.         else if (z/typical <= tolerance)
  237.             equal = true;
  238.     }
  239.  
  240.     return equal;
  241. }
  242.  
  243.  
  244. //---------------------------------------------------------------
  245. //
  246. // SetRandomSeed
  247. //
  248. //---------------------------------------------------------------
  249. void SetRandomSeed(long seed)
  250. {
  251.     sFlipSeed = gb_init_rand(seed);        
  252. }
  253.  
  254.  
  255. //---------------------------------------------------------------
  256. //
  257. // Random (long)
  258. //
  259. //---------------------------------------------------------------
  260. long Random(long max)
  261. {
  262.     return (long) gb_unif_rand(max);                         
  263. }
  264.  
  265.  
  266. //---------------------------------------------------------------
  267. //
  268. // Random (double)
  269. //
  270. //---------------------------------------------------------------
  271. double Random(double max)
  272. {
  273.     long n = gb_unif_rand(LONG_MAX); 
  274.     
  275.     double x = max*scalb(n, -31);
  276.         
  277.     return x;                         
  278. }
  279.  
  280.  
  281. //---------------------------------------------------------------
  282. //
  283. // Random (long, long)
  284. //
  285. //---------------------------------------------------------------
  286. long Random(long min, long max)
  287. {
  288.     long range = max - min;
  289.     
  290.     return min + gb_unif_rand(range);
  291. }
  292.  
  293.  
  294. //---------------------------------------------------------------
  295. //
  296. // Random (double, double)
  297. //
  298. //---------------------------------------------------------------
  299. double Random(double min, double max)
  300. {
  301.     double range = max - min;
  302.     
  303.     return min + Random(range);
  304. }
  305.  
  306. #pragma mark -
  307.  
  308. // ===================================================================================
  309. //    Unit Test
  310. // ===================================================================================
  311.  
  312. //---------------------------------------------------------------
  313. //
  314. // GetRandom
  315. //
  316. // A function to help test the random number testers. It does this
  317. // by providing some low quality random number generators that will
  318. // hopefully do poorly on our tests.
  319. //
  320. //---------------------------------------------------------------
  321. #if DEBUG
  322. static short GetRandom(short max)
  323. {
  324. #if 1
  325.     // Use the real random number generator.
  326.     short x = Random(max);
  327.     
  328. #elif 0
  329.     // Use a real lousy random number generator.
  330.     gRandomSeed = Abs(7*gRandomSeed) % 211;    
  331.     long y = Abs(gRandomSeed) % 100;                    // y is in [0, 99]                
  332.     short x = (short) (max*y/100);                        
  333.  
  334. #elif 0
  335.     // Return the next integer in the range.
  336.     short x = (short) (Abs(gRandomSeed++) % max);
  337. #endif
  338.  
  339.     return x;                         
  340. }
  341. #endif    // DEBUG
  342.  
  343.  
  344. //---------------------------------------------------------------
  345. //
  346. // PassesChiSqr
  347. // 
  348. // See Algorithms by Sedgewick.
  349. //
  350. //---------------------------------------------------------------
  351. #if DEBUG
  352. static bool PassesChiSqr(short max)
  353. {
  354.     long i;
  355.     
  356.     const long kNoIterations = 10L*max;
  357.     long* count = new long[max];
  358.     
  359.     for (i = 0; i < max; i++)
  360.         count[i] = 0;
  361.     
  362.     for (i = 0; i < kNoIterations; i++) {
  363.         short x = GetRandom(max);
  364.         if (x >= 0 && x < max)
  365.             count[x] += 1;
  366.         else
  367.             TRACE("Random(%d) returned %d.¥n", max, x);
  368.     }
  369.     
  370.     long t = 0;
  371.     for (i = 0; i < max; i++) 
  372.         t += count[i]*count[i];
  373.         
  374.     long chiSqr = (max*t/kNoIterations) - kNoIterations;
  375.     long toler = (long) (2*sqrt((float) max));
  376.     
  377. #if 0
  378.     TRACE("Chi-square for Random(%d) = %d¥n", max, chiSqr);
  379.     TRACE("It should be within %d of %d¥n", toler, max);
  380. #endif
  381.     
  382.     delete [] count;
  383.     
  384.     return Abs(chiSqr - max) <= toler;
  385. }
  386. #endif    // DEBUG
  387.  
  388.  
  389. //---------------------------------------------------------------
  390. //
  391. // TestSpread
  392. //
  393. // Here we check to see if the values returned by the random
  394. // number generator are reasonably spread out. We do this by
  395. // using the chi-squared test. This can fail about once in ten
  396. // times, so we'll write an error message if it fails more than
  397. // twice.
  398. //
  399. //---------------------------------------------------------------
  400. #if DEBUG
  401. static void TestSpread(short max)
  402. {
  403.     short failCount = 0;
  404.     
  405.     for (short i = 0; i < 10; i++)
  406.         if (!PassesChiSqr(max))
  407.             failCount++;
  408.             
  409.     if (failCount > 2) {
  410.         TRACE("Random(%d) failed the chi-squared test %d times.¥n", max, failCount);
  411.         TRACE("Note that failing the test once is perfectly normal.¥n¥n");
  412.     }
  413. }
  414. #endif    // DEBUG
  415.  
  416.  
  417. //---------------------------------------------------------------
  418. //
  419. // HasACycle    
  420. //
  421. // Writes a message if cycleLen values starting at startIndex
  422. // form a cycle (ie the next set of cycleLen values is the same
  423. // as the first). 
  424. //
  425. //---------------------------------------------------------------
  426. #if DEBUG
  427. static bool HasACycle(short values[], short noValues, short cycleLen, short startIndex, short max)
  428. {
  429.     short entry = 0;                    // index into cycle
  430.     short repeatCount = 0;                // number of times cycle has been repeated
  431.     short index = (short) (startIndex + cycleLen);// index into next value to be checked
  432.     bool cycling = true;                // turns false when cycle is broken
  433.  
  434.     while (index < noValues && cycling) {
  435.         if (values[startIndex + entry++] != values[index++])
  436.             cycling = false;
  437.         else if (entry == cycleLen) {
  438.             entry = 0;
  439.             repeatCount++;
  440.         }
  441.     }
  442.     
  443.     if (repeatCount > 0) {
  444.         TRACE("A cycle of length %d was repeated %d", cycleLen, repeatCount);
  445.         TRACE(" times with Random(%d).¥n", max);
  446.     }
  447.     
  448.     return repeatCount > 0;
  449. }
  450. #endif    // DEBUG
  451.  
  452.  
  453. //---------------------------------------------------------------
  454. //
  455. // TestForCycles
  456. //
  457. // A good random number generator should generate numbers that
  458. // are independant of the numbers it cranked out before. In this
  459. // test we check to see if the random number generator degenerates
  460. // into cyclic behavior. For example, a bad random number generator
  461. // might fall into a cycle of length 4 for Random(3): 011201120112.
  462. //
  463. // The code below checks for cycles of length max and up. (If the
  464. // random number generator falls into a cycle with length smaller
  465. // than max it will fail the frequency test above). 
  466. //
  467. //---------------------------------------------------------------
  468. #if DEBUG
  469. static void TestForCycles(short max)
  470. {
  471.     const short kNoIterations = 1000;
  472.     short values[kNoIterations];
  473.     
  474.     for (short i = 0; i < kNoIterations; i++)
  475.         values[i] = GetRandom(max);
  476.         
  477.     const short kMaxCycle = kNoIterations/3;
  478.     for (short cycleLen = max; cycleLen < kMaxCycle; cycleLen++)
  479.         for (short startIndex = 0; startIndex < kNoIterations - cycleLen; startIndex++)
  480.             if (HasACycle(values, kNoIterations, cycleLen, startIndex, max))
  481.                 break;
  482. }
  483. #endif    // DEBUG
  484.  
  485.  
  486. //---------------------------------------------------------------
  487. //
  488. // TestRandom
  489. //
  490. //---------------------------------------------------------------
  491. #if DEBUG
  492. static void TestRandom()
  493. {
  494.     while (gb_flip() < 0)                    // Random
  495.         ;
  496.     TestSpread(2);             
  497.     TestSpread(10); 
  498.     TestSpread(105); 
  499.     TestForCycles(10); 
  500.     TestForCycles(23); 
  501.     TestForCycles(100); 
  502.     
  503.     // ･･･ﾊneed to test FlipCoin as well
  504.     
  505.     TRACE("Completed random test.¥n¥n");
  506. }
  507.  
  508. static TUnitTestRegistrar sRandReg("Random", TestRandom);
  509.  
  510. #endif    // DEBUG
  511.